Newton-Cotes integration for approximating Stieltjes (generalized Euler) constants

In the Laurent expansion ζ(s, a) = 1 s− 1 + ∞ ∑ k=0 (−1)γk(a) k! (s− 1) , 0 < a ≤ 1, of the Riemann-Hurwitz zeta function, the coefficients γk(a) are known as Stieltjes, or generalized Euler, constants. [When a = 1, ζ(s, 1) = ζ(s) (the Riemann zeta function), and γk(1) = γk.] We present a new approach to high-precision approximation of γk(a). Plots of our results reveal much structure in the growth of the generalized Euler constants. Our results when 1 ≤ k ≤ 3200 for γk, and when 1 ≤ k ≤ 600 for γk(a) (for a such as 53/100, 1/2, etc.) suggest that published bounds on the growth of the Stieltjes constants can be much improved, and lead to several conjectures. Defining g(k) = sup0 1. We also conjecture that limk→∞ ( γk(1/2) + γk ) /γk = 0, a special case of a more general conjecture relating the values of γk(a) and γk(a+ 1 2 ) for 0 < a ≤ 1 2 . Finally, it is known that γk = limn→∞{ ∑n j=2 log j j − log k+1 n k+1 } for k = 1, 2, . . . . Using this to define γr for all real r > 0, we conjecture that for nonintegral r, γr is precisely (−1)r times the r-th (Weyl) fractional derivative at s = 1 of the entire function ζ(s)− 1/(s− 1)− 1. We also conjecture that g, now defined for all real arguments r > 0, is smooth. Our numerical method uses Newton-Cotes integration formulae for very high-degree interpolating polynomials; it differs in implementation from, but compares in error bounding to, Euler-Maclaurin summation based methods.